The Dynamics of Structural Transitions in Sodium Chloride Clusters 
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In recent experiments on sodium chloride clusters structural transitions between nanocrystals 
with different cuboidal shapes were detected. Here we presents results for the thermodynamics and 
dynamics of one of these clusters, (NaCl)35Cl~. As the time scales for the structural transitions 
can be much longer than those accessible by conventional dynamics simulations, we use a master 
equation to describe the probability flow within a large sample of potential energy minima. We 
characterize the processes contributing to probability flow between the different nanocrystals, and 
obtain rate constants and activation energies for comparison with the experimental values. 



I. INTRODUCTION 

Understanding the relationship between the potential 
energy surface (PES), or energy landscape, and the dy- 
namics of a complex system is a major research effort in 
the chemical physics community. One particular focus 
has been the dynamics of relaxation from an out-of equi- 
librium starting configuration down the energy landscape 
to the state of lowest free energy, which is often also the 
global minimum of the PES.S 

The possible relaxation processes involved can be 
roughly divided into two types. The first is relaxation 
from a high-energy disordered state to a low-energy or- 
dered state, and examples include the folding of a protein 
from the denatured state to the native structure, and the 
formation of a crystal from the liquid. 

The second kind of relaxation process, which is the fo- 
cus of the work here, is relaxation from a low-energy, but 
metastable, state to the most stable state. This second 
situation often arises from the first; the initial relaxation 
from a disordered state can lead to the population of 
a number of low-energy kinetically accessible configura- 
tions. The time scales for this second relaxation process 
can be particularly long because of the large free energy 
barriers that can separate the states. 

Some proteins provide an instance of this second type 
of relaxation. Often, as well as a rapid direct path from 
the denatured state to the native structure, there is also 
a slower path which passes through a low-energy kinetic 
trap.Bl3 As this trapping is a potential problem for pro- 
tein function, the cell has developed its own biochemical 
machinery to circumvent it. For example, it has been 
suggested that the chaperonin, GroEL, aids protein fold- 
ing by unfolding those protein molecules which get stuck 
in a trapped stateD 

There are also a growing number of examples of 
this second type of relaxation involving clusters. For 
Lennard-Jones (LJ) clusters there are a number ofj-sizes 
for which the global minimum is non-icosahedral.tj For 
example, for LJ38 the global minimum is a face-centred- 
cubic truncated octahedron, but relaxation down the 
PES almost always leads to a low-energy icosahedral min- 
imum. Similarly, for LJ75 the icosahedral states act as 
a trap preventing relaxation to the Marks decahedralcl 



global minimum. A similar competition between face- 
centred-cubic or decahedral and icosahedral .structures 
has recently been observed for metal clusters.! 

For these clusters unbiased global optimization is dif- 
ficult because the icosahedral states act as an effective 
trap. More generally, kinetic traps are one of the major 
problems for a global optimization algorithm. Therefore, 
much research has focussed on decreasing the 'life-time' 
of such traps. For example, some methods simulate a 
non-Boltzmann ensemble that involves increasetLfluctu- 
ations, thus making barrier crossing more likelylm Other 
algorithms transform the energy landscape in a way that 
increases the temperature range where the global mini- 
mum is populated, thus allowing one to choose conditions 
where the free /energy barriers relative to the thermal en- 
ergy are lower 

Recently, clear examples of trapping associated with 
structural transitions have emerged from experiments on 
NaCl clusters. These clusters have only one energeti- 
cally favourable morphology: the magic numbers that 
appear in mass spectra correspond to .cubaidal frag- 
ments of the bulk crystal (rocksalt) lattice,ll3~EJ hence the 
term nanocrystals. Indirect structural information comes 
from the experiments of Jarrold and coworkers which 
probe the mobility of size-selected cluster ions. For most 
(NaCl)jvCl - with N > 30, multiple isomers were de- 
tected which wece-assigned as nanocrystals with different 
cuboidal shapes .Ej The populations in the different iso- 
mers were not initially equilibrated, but slowly evolved, 
allowing rates and activation energies for the structural 
transitions between the nanocrystals to be obtained.li3 

In a previous paper we identified the mechanisms of 
these structural transitions by extensively searching the 
low-energy regions of the PES of one of these clus- 
ters, (NaCl^Cl - , in order to obtain paths linking the 
different cuboidal morphologies. ell The key process in 
these transitions is a highly cooperative rearrangement 
in which two parts of the nanocrystal slip past one an- 
other on a {110} plane in a (110) direction. 

Here we continue our examination of the struc- 
tural transitions by investigating the dynamics of 
(NaCl)35Cl~. Given the long time scales for which the 
clusters reside in metastable forms, it is not feasible to 
probe the transitions with conventional dynamics simu- 
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lations. Instead, we use a master equation that describes 
the probability flow between the minima on the PES. 
This method has the advantage that we can relate the 
dynamics to the topography of the PES . In this paper 
we are particularly concerned with obtaining activation 
energies for the structural transitions: firstly, in order to 
compare with experiment, and secondly, to understand 
how the activation energy for a process that involves a 
series of rearrangements and a large number of possible 
paths is related to the features of the energy landscape. 
In section |l| we outline our methods, and then in Section 



III 



after a brief examination of the topography of the 
PES and the thermodynamics, we present our results for 
the dynamics of the structural transitions. 



(NaCl) 35 Cl-.0 This sampling was performed by repeat- 
edly stepping across the PES from minimum to minimum 
via transition states, thus giving a connected set of sta- 
tionary points. We biased this search to probe the low- 
energy regions of the PES either by using a Metropo- 
lis criterionEll|-to|-decide whether to accept a step to a 
new minimumE3~E-3 or by systematically performing tran- 
sition stata-searches from the lower energy minima in 
the sampleE§E3 Thus, although our samples of station- 
ary points only constitute a tiny fraction of the total 
number on the PES, we have a good representation of 
the low-energy regions that are relevant to the structural 
transitions of (NaCl)35Cl~. 



II. METHODS 



A. Potential 



The potential that we use to describe the sodium chlo- 
ride clusters is the Tosi-Fumi parameterization of the 
Coulomb plus Born-Mayer formtS: 



Aije~ 



i/p 



where qt is the charge on ion i, is the distaace between 
ions i and j and Aij and p are parameters^. Although 
simple, this potential provides a reasonable description 
of the interactions. For example, in a previous study we 
compared the global minima of (NaCl)]vCl~ (A < 35) 
for this potential with those for a more complex poten- 
tial derived by Welch, et al. which also includes terms 
due to polarization.EJa Most p£ the global minima were 
the same for both potentialsO Given some of the other 
approximations we use in this study, the small advantages 
gained by using the Welch potential do not warrant the 
considerable additional computational expense. 

We should also note that a well-known problem as- 
sociated with the above family of potentials for the al- 
kali halides is that they never predict the CsCl structure 
to be the most stable. This problem arises because the 
potentials do not allow the properties of pi ion to be 
dependent on the local ionic environment E3 This defi- 
ciency should not greatly affect the relative energies of 
the low-lying (NaCl^Cl - minima because they all have 
the same rock-salt structure, but it may affect the bar- 
riers for rearrangements where some ions experience a 
different local environment at the transition state. 



B. Searching the potential energy surface 

The samples of 3518 minima and 4893 transition 
states that we use here were obtained in our previous 
study on the mechanisms of the structural transitions for 



C. Thermodynamics 

We used two methods to probe the thermodynamics: 
first, conventional Monte Carlo simulations and second, 
the superposition method. The latter is a technique to 
obtain thq thermodynamics of a system from a sample of 
minimaB3 It is based on the idea that all of configura- 
tion space can be divided upifito the basins of attraction 
surrounding each minimumfia The density of states or 
partition function can then be written as a sum over all 
the minima on the PES, e.g. Z — Zi, where Z^ is the 
partition function of minimum i. 

The limitations of the superposition method are that 
the Zi are not known exactly and that, for all but the 
smallest systems, the total number of minima on the PES 
is too large for us to characterize them all. However, the 
harmonic expression for Zi leads to a reasonable descrip- 
tion of the thermodynamics.^!! Furthermore, anharmonic 
forms are available which allow the thermodynamics of 
larger clusters to be reproduced accurately.^ 

The incompleteness of the sample can be overcome by 
weighting the contributions from the minima in a repre- 
sentative sample.EJ However, this approach is not neces- 
sary in the present study since we are interested in low 
temperature behaviour where the number of thermody- 
namically relevant minima is still relatively small. Fur- 
thermore, in this temperature regime the superposition 
method has the advantage that it is unaffected by large 
free energy barriers between low-energy minima which 
can hinder the determination of equilibrium thermody- 
namic properties by conventional simulation. 

Here we use the harmonic form of the superposition 
method, because we later use the harmonic approxima- 
tion to derive rate constants (reliable anharmonic expres- 
sions for the rate constants are not so readily available). 
The partition function is then 



E 



-0Ei 



(i) 



where f3 = 1/kT, Ei is the energy of minimum i, is the 
geometric mean vibrational frequency of i, k — 3 A — 6 
is the number of vibrational degrees of freedom and rii 
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is the number of permutational isomers of i. rii is given 
by 2N\/hi, where hi is the order of the point group of i. 
From this expression thermodynamic quantities such as 
the heat capacity, C v , can be obtained by the application 
of standard thermodynamic formulaeEH The superposi- 
tion method also allows us to examine the contributions 
of particular regions of configuration space to the the 
thermodynamics. For example, the probability that the 
system is in region A is given by 



PA 



Zi 
Z' 



(2) 



where the sum runs over those minima that are part of 
region A. 

D. Dynamics 

The master equation approach^ is increasingly be- 
ing used to describe the inter-minimum dynamics on a 
multi-dimensprLal PES with, applications, ta, for exam- 
ple, clusters,tirE-3 glasses,E§E3 proteinsEZn^l and ideal- 
ized model landscapes.HI The master equation is defined 
in terms of P(t) = {Pi(t)}, the vector whose components 
are the ensemble-average probabilities that the system is 
associated with each minimum at time t: 



dPjjt) 
dt 



^[kijPjW-kjiPiit)], 



(3) 



where kij is the rate constant for transitions from mini- 
mum j to minimum i. Defining the matrix 



the zero eigenvalue, all Xj are negative. Therefore, as 
t — > oo the contribution of these modes decays exponen- 
tially to zero and P — > V eq . 

To apply Equation (|s|) we must first diagonalizc W. 
The computer time required for this procedure scales as 
the cube of the size of the matrix and the memory re- 
quirements scale as the square. Therefore, it is advanta- 
geous for the matrix W to be as small as possible. For 
this reason we recursively removed those minima that are 
only connected to one other minimum; these 'dead-end' 
minima do not contribute directly to the probability flow 
between different regions of the PES. After pruning our 
samples have 1624 minima and 2639 transition states. To 
test the effect of this pruning, we performed some calcu- 
lations using both the full and the pruned samples. The 
effect on the dynamics of the structural transitions was 
negligible. 

As the temperature of a system is decreased the spread 
of eigenvalues can increase rapidly. When the ratio of the 
largest to smallest eigenvalues is of the order of the ma- 
chine precision of the computer, the accuracy of the ex- 
treme eigenvalues can become degraded by rounding er- 
rors. We encountered these problems below 275K. With- 
out pruning the samples these numerical difficulties are 
more pronouncedfia 

We model the rate constants, which arer^eeded as in- 
put to Equation (^), using RRKM theoryc3 in the har- 
monic approximation. Therefore, the rate constant for a 
transition from minimum i to minimum j via a particular 
transition state (denoted by f ) is given by 



4(T) = -^- T |^ y exp(-(4-^)/fcT). 



(6) 



(4) 



III. RESULTS 



allows Equation (^|) to be written in matrix form: 
dP(t)/dt = WP(t). 

If the transition matrix W cannot be decomposed into 
block form, it has a single zero eigenvalue whose corre- 
sponding eigenvector is the equilibrium probability dis- 
tribution, V eq . As a physically reasonable definition for 
the rate constants must obey detailed balance at equilib- 



rium, i.e. WijPp 



WaP° q , the solution of the master 



equation can be expanded in terms of a complete set of 
eigenf unctions of the symmetric matrix, W, defined by 
Wij = (P^/P^^Wij. The solution is 



3=1 




(5) 



where u\ is component i of the j eigenvector of W 
and Xj is the j th eigenvalue. 

The eigenvalues of W and W are identical and the 
eigenvectors are related by = y ' P^ q ' . Except for 



A. Topography of the (NaCl) 35 Cr PES 

In our previous study of (NaCl^Cl - we found that 
the low-energy minima all had rock-salt structures. The 
different minima have four basic shapes: an incomplete 
5x5x3 cuboid, a 6 x 4 x 3 cuboid with a single vacancy, a 
8x3x3 cuboid with a single vacancy and an incomplete 
5x4x4 cuboid. The lowest-energy minimum for each 
of these forms is shown in Figure [l| 

In the experiments on (NaCl^Cl^ by Jarrold and 
coworkers the three peaks that were resolved in the ar- 
rival time distribution were assigned on the basis of cal- 
culated mobilities as 5x5x3, 5x4x4 and 8x3x3 
nanocrystalsJiaij However, when the 6x4x3 nanocrys- 
tal is also considered, better agreement between the cal- 
culated and observed mobilities can be obtained by as- 
signing the three experimental peaks to the 6x4x3, 
5x5x3 ap^n-Si x 3 x 3 nanocrystals in order of increasing 
drift time£3EJ This reassignment is also in better agree- 
ment with the thermodynamics since the clusters convert 
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to (what is now assigned as) the 5x5x3 nanocrystal 
as time progresses, Ej indicating that this structure has 
the lowest free energy. In our calculations a 5 x 5 x 3 
isomer is the global potential energy minimum, and the 
5x5x3 nanocrystal is al ways m ore stable than the other 
nanocrystals (See Section III B ). 
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FIG. 1. The lowest-energy 5 x 5 x 3 (A), 6 x 4 x 3 (D) and 
8 x 3 x 3 (I) nanocrystals and the two lowest-energy 5x4x4 
nanocrystals (L and O). The sodium ions are represented 
by the smaller, darker circles and the chloride ions by the 
larger, lighter circles. The letter gives the energetic rank of 
the minimum when labelled alphabetically. 

Disconnectivity graphs provide a way of visualiz- 
ing an energy landscape that is particularly useful for 
obtaining insight into dynamics, and have prcyiously 
been applied te> a number of protein modelsEZlc3~EJ and 
clusters. c3e3 _ c3 The graphs are constructed by perform- 
ing a 'superbasin' analysis at a series of energies. This 
analysis involves grouping minima into disjoint sets, 
called superbasins, whose members are connected by 
pathways that never exceed the specified energy. At each 
energy a superbasin is represented by a node, and lines 
join nodes in one level to their daughter nodes in the 
level below. Every line terminates at a local minimum. 
The graphs therefore present a visual representation of 
the hierarchy of barriers between minima. 

The disconnectivity graph for (NaCl^sCl - is shown in 
Figure ||. The barriers between minima with the same 
cuboidal form are generally lower than those between 
minima that have different shapes. Therefore, the dis- 
connectivity graph splits into funnels corresponding to 
each cuboidal morphology. (A funnel is a set of downhill 
pathways that converge on a single low-energy muiimum 
or a set of closely-related low-energy minima.cilEj In a 
disconnectivity graph an ideal funnel is represented by 
a single tall stem with lines branching directly from it, 
indicating the progressive exclusion of minima as the en- 
ergy is decreased. c3) The separation is least clear for the 
5x4x4 minima because of the large number of different 



ways that the nine vacant sites can be arranged. For ex- 
ample, these vacancies are organized very differently in 
the two lowest-energy 5x4x4 minima (Figure [l]) , and in 
fact the barrier between minimum O and the low-energy 
6x4x3 isomers is lower than the barrier between O 
and minimum L. Therefore, the minima associated with 
O form a sub-funnel that splits off from the 6x4x3 
funnel, rather than being directly connected to the main 
5x4x4 funnel. 

The disconnectivity graph shows that the barriers be- 
tween the 5x5x3,6x4x3 and 5x4x4 nanocrystals 
are of similar magnitude, while the 8x3x3 minima are 
separated from the rest by a considerably larger barrier. 
The values of some of the barrier heights are given in 
Table |. 



-264.2 ,£/eV 
-264.3 
-264.4 
-264.5 
-264.6 
-264.7 
-264.8 
-264.9 
-265.0 
-265.1 
-265.2 
-265.3 
-265.4 
-265.5 
-265.6 
-265.7 
-265.8 

FIG. 2. Disconnectivity graph for (NaCl) 35 Cl". Only 
branches leading to the lowest 200 minima are shown. The 
branches for the twenty lowest-energy minima are labelled 
alphabetically in order of their energy. Some of the funnels 
and sub-funnels have also been labelled by their associated 
cuboidal shape. 

The disconnectivity graph is also helpful for in- 
terpreting thft— 1 (NaCl)35Cl~ dynamics observed in 
experiments.ll3 : t3 In the formation process it is likely that 
a high-energy configuration is initially generated. The 
cluster then relaxes to one of the low-energy nanocrys- 
tals. Simulations for potassium chloride clusters indicate 




8x3x3 
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that this relaxation is particularly rapid for alkali halides 
because of the low barriers (relative to the energy,-dit| 
ference between the minima) for downhill pathways. lJ'EHI 
However, there is a separation of time scales between this 
initial relaxation and the conversion of the metastable 
nanocrystals to the one with the lowest free energy, the 
large barriers between the different cuboids make them 
efficient traps. 




300 400 500 600 700 800 900 

temperature / K 



FIG. 3. The temperature dependence of (a) the total en- 
ergy and (b) the radius of gyration, R g , for four series of MC 
runs starting in the lowest-energy minima of the 5x5x3 
(diamonds), 6x4x3 (plus signs), 8x3x3 (squares) and 
5x4x4 (crosses) nanocrystals. Each point is the average 
value in a 10 6 cycle Monte Carlo run, where each run was 
initiated from the final geometry in the previous lower tem- 
perature run. The error bars in (b) represent the standard 
deviation of the R g probability distributions. 



B. Thermodynamics of (NaCl) 36 Cl~ 

Some thermodynamic properties of (NaCl^sCl - are 
shown in Figures || and ||. The caloric curve shows a fea- 
ture at ^700K which indicates melting (Figure |a); the 
melting temperature is depressed relative to the bulk due 



to the cluster's finite size. The effects of the barriers be- 
tween nanocrystals are apparent in the MC simulations. 
The radius of gyration, R g , provides a means of differen- 
tiating the nanocrystals. ft can be seen from the plot of 
R g for simulations started in the lowest-energy minima 
of the four cuboidal forms that each simulation is stuck 
in the starting structure (Figure [|b) up to temperatures 
close to melting, implying that there are large free energy 
barriers between the nanocrystals. 



(a) 




206 -I 1 1 1 1 

200 400 600 800 1000 

temperature / K 

FIG. 4. Thermodynamic properties of (NaCl^sCP com- 
puted using the harmonic superposition method with the 
pruned sample of minima, (a) Equilibrium occupation prob- 
abilities of the nanocrystals, as labelled, (b) Heat capacity. 

These free energy barriers prevent an easy determina- 
tion of the relative stabilities of the different nanocrystals 
by conventional simulations. Therefore, we use the super- 
position method to examine this question. First we as- 
sign the fifty lowest-energy minima to one of the cuboidal 
forms by visual inspection of each structure. Then, using 
these sets as definitions of the nanocrystals in Equation 
, we calculate the equilibrium probabilities of the clus- 
ter being in the different cuboidal morphologies as a func- 
tion of temperature. It can be seen from Figure |^a that 
the 5x5x3 nanocrystal is most stable up until melting. 
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The 6x4x3 nanocrystal also has a significant probabil- 
ity of being occupied. However, the probabilities for the 
8x3x3 and 5x4x4 nanocrystals are always small. The 
onset of the melting transition is indicated by the rise of 
p rcs t in Figure and by the peak in the heat capacity 
(Figure |Jb). However, this transition is much too broad 
and the heat capacity peak occurs at too high a tempera- 
ture because the incompleteness of our sample of minima 
leads to an underestimation of the partition function for 
the liquid-like minima. Given these expected failings of 
the superposition method at high temperature when the 
partition functions of the minima are not reweighted, we 
restrict our dynamics calculations to temperatures below 
600K. 

Although the probabilities of being in the different 
morphologies show little variation at low temperature, 
there are significant changes in the occupation probabil- 
ities of specific minima. For example, the small low tem- 
perature peak in the heat capacity is a result of a redis- 
tribution of probability amongst the low-energy 5x5x3 
minima; the third lowest-energy minimum becomes most 



populated. It is also interesting to note that the second 
lowest-energy 5x4x4 minimum (O) becomes more stable 
than the lowest-energy 5x4x4 minimum (L) for tem- 
peratures above approximately 220K. Both these changes 
are driven by differences in vibrational entropy. 

C. Dynamics of (NaCl) 35 Cr 

Some examples of the interfunnel dynamics that we 
find on solution of the master equation are depicted in 
Figure ||. The time scales involved are much longer than 
those accessible by conventional simulations. 

The dynamics of relaxation to equilibrium depend 
significantly on the starting configuration. When the 
lowest-energy 6x4x3 minimum is the initial configura- 
tion there is a small transient population in 5 x 4 x 4 min- 
ima before the system adopts a 5 x 5 x 3 structure. This 
is consistent with the lowest-energy pathway that was 
found between the two nanocrystalsp-jt passes through 
some intermediate 5x4x4 minima.ED 




10'' 



time / s 



0.0001 0.001 
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FIG. 5. Time evolution of the occupation probabilities of the different nanocrystals at T=400K when the cluster is 
initially in minimum (a) D, (b) I, (c) L and (d) O. 
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The relaxation to equilibrium is much slower when the 
cluster starts from the lowest-energy 8x3x3 minimum. 
This is a result of the large barrier to escape from this 
funnel (Figure ||) . The probability flow out of the 8x3x3 
funnel leads to a simultaneous rise in the occupation 
probabilities of both the 5x5x3 and 6x4x3 nanocrys- 
tals towards their equilibrium values, even though the 
lowest-energy pathway out of the 8x3x3 funnel di- 
rectly connects it to the low-energy 6x4x3 minima. 
This occurs because the time scale for interconversion of 
these latter two nanocrystals is much shorter than that 
for escape from the 8x3x3 funnel. However, we do find 
evidence that the cluster first passes through the 6x4x3 
minima if we examine the probabilities on a log-scale. At 
times shorter than that required for local equilibrium be- 
tween the 5x5x3 and 6x4x3 minima the occupation 
probability for the 6x4x3 minima is larger. 

As the two lowest-energy 5x4x4 minima are well- 
separated in configuration space we considered relaxation 
from both these minima. In both cases, there is a large 
probability flow into the 6x4x3 minima, which is then 
transferred to the 5x5x3 funnel on the same time scale 
as when initiated in the 6x4x3 funnel. However, the 
time scale for the build-up of population in the 6x4x3 
minima depends on the initial configuration. Probability 
flows more rapidly and directly into the 6x4x3 min- 
ima when initiated from minimum O, reflecting the low 
barriers between these minima (Figure ^|). For relaxation 
from minimum L there are two active pathways, leading 
to an increase in the population of both the 5x5x3 and 
the 6x4x3 minima. The direct path into the 5x5x3 
funnel has the lower barrier (Table ||) but is long (96. 7A), 
and so has a smaller rate than the path into the 6x4x3 
funnel, which has a slightly higher barrier (by 0.05 eV). 
The small shoulder in the occupation probability of the 
5x5x3 minima occurs in the time range when the occu- 
pation probability of the 5x4x4 minima has reached a 
value close to zero (thus reducing the contribution from 
the direct path) and when the probability flow out of the 
5x4x4 minima is only just beginning. 

The combination of our thermodynamics and dynam- 
ics results for (NaCl^Cl - enable us to explain why a 
peak associated with the 5x4x4 cuboids was not ob- 
served experimentally. The 5x4x4 minima have a 
shorter lifetime than the other cuboidal forms and have 
a low equilibrium occupation probability. 

Another way to analyse the dynamics is to exam- 
ine how local equilibration progresses towards the point 
where global equilibrium has been obtained. To accom- 
plish this we define two minima to be in local equilibrium 
at the time when 



cq I 



< e 



(7) 



is obeyed for all later times. In the present work we set 
e = 0.01, i.e. the two minima are within 1% of equilib- 
rium. Using this definition we can construct equilibra- 



tion graphs, in which nodes occur when two (groups of) 
minima come into local equilibrium]^ 

We show an example of an equilibration graph in Fig- 
ure [| Equilibration first occurs at the bottom of the 
6x4x3 and 5x5x3 funnels between those minima 
that are connected by low barrier paths, then progresses 
to minima with the same cuboidal shape but which are 
separated by larger barriers, and finally occurs between 
the funnels. The order of intcrfunncl equilibrium agrees 
with the time scales that we observe in the time evolu- 
tion of the occupation probabilities of the nanocrystals 
(Figure ^|). Minimum O, then minimum L, come into 
equilibrium with the 6x4x3 funnel. Then, the 5x5x3 
and 6x4x3 funnels reach local equilibrium. Finally, the 
8x3x3 funnel reaches equilibrium with the rest of the 
PES. As one of the major determinants of the time scale 
required for local equilibrium is the height of the barriers 
between minima, it is unsurprising that the equilibration 
graph reflects the structure of the disconnectivity graph 
(Figure |). 
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FIG. 6. Equlibration graph showing how the system pro- 
gresses towards equilibrium at T=400K when the cluster is 
initially in minimum D. The lines join when minima come 
into equilibrium with each other. The vertical scale and the 
horizontal position of the ends of the lines are chosen for clar- 
ity. The ends are labelled by the letter for the corresponding 
minimum. 

In the experiments on (NaCl^Cl - rate constants and 
activation energies were obtained for the conversion of 
the 6 x 4 x 3-and 8x3x3 nanocrystals into the 5x5x3 
nanocrystalilS It would, therefore, be useful if we could 
extract rate constants for the different interfunnel pro- 
cesses from the master equation dynamics. 

For a two-state system, where A # B and k+ and fc_ 
are forward and reverse rate constants, respectively, it 
can be shown that 



In 



p<=q 



p a (o)-p c ; 



= -(fc+ + k-)t, 



(8) 
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and the equivalent expression for B are obeyed. M This 
is a standard result for a first-order reaction. This ex- 
pression will also hold for the rate of passage between 
two funnels in our multi-state system if the interfunncl 
dynamics are the only processes affecting the occupation 
probabilities of the relevant funnels, and if the interfun- 
nel dynamics cause the occupation probabilities of the 
two funnels to converge to their equilibrium values. 

In Figure fj] we test the above expression by apply- 
ing it to the interconversion of 6 x 4 x 3 and 5x5x3 
nanocrystals. The two lines in the graph converge to 
the same plateau value, before both falling off beyond 
0.001s. This plateau corresponds to the time range for 
which the interfunnel passage dominates the evolution of 
the probabilities for the two funnels. At shorter times, 
when the occupation probabilities for the two funnels are 
still close to their initial values, there are many other 
contributing processes. At longer times the probabilities 
are both very close to their equilibrium values, and the 
slower equilibration with the 8x3x3 funnel dominates 
the probability evolution. From the plateau in Figure 
we obtain k + + k" = 5320 s _1 . The individual rate 
constants can be obtained by using the detailed balance 



relation: fc + P^ q 



6x4x3 



5x5x3, 



I0" 7 10"" 10 J 0.0001 

time / s 

FIG. 7. By following (1/t) log[(pi(t) - pT)KpM - pT)\ 
at T=400K when minimum D is the initial configuration, 
this graph tests the applicability of Equation (^) to the in- 
terconversion of the 5x5x3 and 6x4x3 nanocrystals. 

The application of Equation (||) to escape from the 
8x3x3 funnel also leads to a range of t where there 
is a well-defined plateau. However, this approach works 
less well for the interconversion of 5 x 4 x 4 and 6x4x3 
minima (Figures ||a and b) . This is because the assump- 
tion that no other processes contribute to the probabil- 
ity evolution of the two funnels is obeyed less well, and 
because the occupation probabilities do not converge to 
their equilibrium values, but near to the equilibrium val- 
ues that would obtain if the 5x5x3 minima were ex- 
cluded. Nevertheless, approximate values of k + and k" 
can be obtained. 



Diagonalization of the matrix W produces a set of 
eigenvalues that give the time scales for a set of char- 
acteristic probability flows. The dynamical processes to 
which the eigenvalues correspond can be identified by 
examining the eigenvectors. Flow occurs between those 
minima for which the corresponding components of the 
eigenvector have opposite sign (Equation (§)). This ob- 
servation forms the basis for Kunz and Berry's net-flow 
index which quantifies the contribution of an eigenvector 
i to flow out of a funnel AOS The index is defined by 



(9) 



The index allows the interfunnel modes to be identified; 
the values of f A and f B for these modes will be large 
and of opposite sign. For example, at T=400K the mode 
with the most 5x5x3 — > 6x4x3 character has 
^5x5x3 = _ .339, / 6x4x3 = 0.331 and A = 5275 s- 1 . 
The eigenvalue is in good agreement with the sum of in- 
terfunnel rates obtained using Equation (^|). 

The extraction of the interfunnel rate in this manner 
is hindered by the fact that the eigenvalues of W can- 
not cross as a function of temperature. Instead, there are 
avoided crossings and mixing of modes. For example, the 
small difference between the two values for the sum of the 
5x5x3^6x4x3 interfunnel rate constants that we 
obtained above is probably due to mixing. The eigenvec- 
tor with which the interfunnel mode mixes gains some 
5x5x3^6x4x3 character; it has / 5x5x3 = -0.014, 
^6x4x3 = 0.017 and A = 6388s" 1 . 



L0 



0.0015 



0.002 



0.0025 0.003 
1 IT 



0.0035 



0.004 



FIG. 8. Arrhenius plots for the rates of interfunnel 
passage. The data points are derived from the eigenval- 
ues of the matrix W and the lines are fits to the form 
k — AexTp(~E a /kT). In ascending order in the figure the 
lines are for the 5x5x3^8x3x3, 8x3x3^5x5x3, 
5x5x3 ^ 6x4x3, 6x4x3^5x5x3, 6x4x3 -> 5x4x4(0), 
6x4x3-»5x4x 4(L), 5 x 4 x 4(0) -» 6 x 4 x 3 and 
5 x 4 x 4(i) -^6x4x3 processes. 
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By calculating and diagonalizing W at a series of tem- 
peratures we can examine the temperature dependence of 
the interfunnel rate constants. For most of the interfun- 
nel processes, the rate constants that we obtain fit well 
to the Arrhenius form, k = Acxp(—E a /kT), where E a is 
the activation energy and A is the prefactor (Figure ||). 

In our previous study of the interconvcrsion mecha- 
nisms in (NaCl^Cl - we estimated the activation en- 
ergies in order to compare with the experimental val- 
ues. Using the analogy to a simple one-step reaction, 
we equated the activation energy with the difference in 
energy between the highest-energy transition state on 
the lowest-energy path between the relevant nanocrys- 
tals and the, lowest-energy minimum of the starting 
nanocrystalO However, it was not clear how well this 
analogy would work. In Table | we compare these es- 
timates with the activation energies obtained from our 
master equation results. There is good agreement, con- 
firming the utility of the approximate approach. Similar 
agreement has also been found for the interfunnel dy- 
namics of a 38-atom Lennard-Jones cluster .o 

A simple explanation for this correspondence can be 
given. We first label the minima along the lowest-barrier 
path between the two funnels in ascending order and de- 
fine I so that the highest-energy transition state on this 
path lies between minima I — 1 and I. If the minima be- 
hind the highest-energy transition state (i.e. 1 to I — 1) 
are in local equilibrium, then the occupation probability 
of the minimum Z — 1 is given by 

V — oc exp {-{Ei-! - E x )/kT) (10) 
Pi 

Then, if the rate of interfunnel flow, k + pA, is equated to 
the rate of passage between minima I — 1 and I, 

k + pA oc pi-i exp{-{Ej_ ll - Ei-{)/kT) 

c< Pl exp(- (i^ -£i)/fcT), (11) 

Therefore, if the occupation probability for funnel A is 
dominated by the occupation probability for the lowest- 
energy minimum in the funnel, i.e. pa ~ Pi for all T, 
then the activation energy is equal to the energy differ- 
ence between the highest-energy transition state on the 
lowest-barrier path and the energy of the lowest-energy 
minimum in the starting funnel. 

We should note that in the above derivation the inter- 
funnel probability flow is assumed to all pass through a 
single transition state. However, if there is competition 
between two paths, one with a low barrier and a small 
prefactor and one with a larger barrier and a large pref- 
actor, we expect that the low-barrier path would dom- 
inate at low temperature and the high-barrier path at 
high temperature. This behaviour would give rise to a 
interfunnel rate constant with a positive curvature in an 
Arrhenius plot. However, the lines are either straight or 
have a small amount of negative curvature. 

The lack of positive curvature, and the agreement be- 
tween the estimated and the observed activation energies, 



probably indicates that the interfunnel probability flow 
is dominated by paths which pass through the highest- 
energy transition state on the lowest barrier path. It is 
interesting that, on a PES with so many minima and 
transition states, a single transition state can have such 
a large influence on the dynamics. 

At low enough temperature, d(pi / 'p a) / 'dT ~ 0, and so 
the interfunnel barrier height can be measured with re- 
spect to the lowest-energy minimum in the starting fun- 
nel (Equation (|ll"|)). However, as the occupation proba- 
bilities of other minima in the funnel becomes significant 
relative to that for the lowest-energy minimum in the fun- 
nel, the ratio Pi/pa decreases, thus giving the lines in the 
Arrhenius plot their slight negative curvature. In other 
words, the apparent activation energy decreases with in- 
creasing temperature, because the barrier height should 
be measured with respect to some kind of average min- 
imum energy for the funnel, perhaps Ea = J2i£APi^i- 
The negative curvature is most pronounced when the oc- 
cupation probabilities in a funnel change considerably. 
For example, the 5 x 4 x 4{L) -^6x4x3 rate constant 
has the most curvature because minimum L has a par- 
ticularly low vibrational entropy leading to population of 
other minima within that funnel (Figure ^) . 

In Table | we also compare our activation energies to 
the two experimental values. Our values are too large by 
0.24 and 0.49 eV. There are a number of possible sources 
of error. Firstly, the samples of minima and transition 
state provide only an incomplete characterization of the 
PES and so it is possible that the nanocrystals are con- 
nected by undiscovered lower-barrier paths. However, we 
believe that our sample of minima is a good representa- 
tion of the low-energy regions of the PES and consider 
it improbable that undiscovered pathways could account 
for all of the discrepancy. Secondly, in our calculations 
the input rate constants fcj were calculated on the basis 
of the harmonic approximation. Although this is likely 
to have a significant effect on the absolute values of the 
interfunnel rate constants and prefactors (the latter are 
too large compared to the experiment valuestj) , it should 
not have such a significant effect on the activation ener- 
gies. Instead, we consider the most likely source of the 
discrepancy between theory and experiment to be inac- 
curacies in the potential. When polarization is included 
by using the Welch potential, the estimated barriers be- 
come closer to the experimental values (the discrepancies 
are then_0.16 and 0.34 eV) but significant differences still 
remain.EZI For better agreement we may need a potential 
that allows the properties of the ions to depend on the 
local environment. Unfortunately, although such pateru 
tials have been developed for a number of systems ,E-lE3 
one does not yet exist for sodium chloride. 
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IV. CONCLUSION 

The (NaCl^Cl - potential energy surface has a 
multiple-funnel topography. Structurally, the different 
funnels correspond to rocksalt nanocrystals with differ- 
ent cuboidal forms. The large potential energy barriers 
between the funnels causes the time scales for escape from 
metastable nanocrystals to be far longer than those acces- 
sible by conventional dynamics simulations. Therefore, 
we examined the interfunnel dynamics by applying the 
master equation approach to a database of (NaCl)35Cl~ 
minima and transition states. The slowest rate constant 



<<5x 5x3^8x3x3 



= 1.46 x ifr 8 s- 



at 



we obtained was 
T = 275K. 

Using a net flow index we were able to identify the 
eigenvalues of the transition matrix, W, which corre- 
spond to interfunnel probability flow. Thus, we were 
able to obtain rate constants and activation energies for 
the interconversion of the different nanocrystals. One 
particularly interesting finding is that the activation en- 
ergies correspond fairly closely to the potential energy 
differences between the highest-energy transition state 
on the lowest-energy path between two nanocrystals and 
the lowest-energy minimum of the starting nanocrystal. 
This is the result one might expect by a simple extrapo- 
lation from the dynamics of a simple molecular reaction. 
However, it holds despite the multi-step, and potentially 
multi-path, nature of the interfunnel dynamics. The 
question of whether this result is generally true for inter- 
funnel dynamics involving large potential energy barriers 
or reflects some of the particulars of the (NaCl^Cl - sys- 
tem is an interesting subject for further research. We al- 
ready know that this simplification holds for the 38-atom 
Lennard-Jones cluster Ha 
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TABLE I. Activation energies, E a , and prefactors, A, 
for the interfunnel rate constants obtained from the fits in 
Figure |8|. The activation energies are compared to AE, the 
difference in energy between the highest-energy transition 
state on the lowest-energy path directly between the two 
nanocrystals and-ithe lowest-energy minimum of the start- 
ing nanoqfYstal,li3 and experimental values for some of the 
processes." As probability flow out of minimum L leads to 
population of both the 6x4x3 and 5x5x3 nanocrystals 
we give AE for both pathways. 









barrier/eV 


A/s- 1 


From 


To 


E a 


AE 


Expt. 


6x4x3 


5x5x3 


0.770 


0.776 


0.53 ±0.05 


1.43 x 10 13 


5x5x3 


6x4x3 


0.846 


0.840 




2.26 x 10 13 


8x3x3 


5x5x3 


1.055 


1.055 


0.57 ±0.05 


1.05 x 10 15 


5x5x3 


8x3x3 


1.220 


1.211 




3.77 x 10 14 


5 x 4 x 4 L 


6x4x3 


0.651 


0.668 




8.46 x 10 13 


5 x 4 x 4 L 


5x5x3 




0.618 






6x4x3 


5 x 4 x 4 L 


0.822 


0.810 




3.95 x 10 14 


5x5x3 


5 x 4 x 4 L 




0.823 






5x4x406x4x3 


0.568 


0.560 




4.97 x 10 13 


6x4x3 


5 x 4 x 5 O 0.738 0.738 




2.32 x 10 14 
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